Basic usage
from langesim import Simulator
Define the stifness k(t) and center c(t) of the harmonic potential
def k(t):
"""Stifness of the potential.
Varies linearly from ki=0.5 to kf=1.0 in tf=5.0 time
units, then remain constant at kf
"""
tf = 5.0
ki = 0.5
kf = 1.0
if t <= 0.0:
return ki
if t < tf:
return ki + (kf - ki) * t / tf
if t >= tf:
return kf
def c(t):
"""Center of the harmonic potential.
Varies linearly from ci = -1.0 to cf = 1.0 in tf time units, then remains constant.
"""
tf = 5.0
ci = -1.0
cf = 1.0
if t <= 0.0:
return ci
if t < tf:
return ci + (cf - ci) * t / tf
if t >= tf:
return cf
Create the simulator
simulator = Simulator(
tot_sims = 100_000,
dt = 0.001,
tot_steps = 10_000,
snapshot_step = 100,
k = k,
center = c
)
The simulator has not performed any simulations yet:
print(simulator.simulations_performed)
0
Run a simulation with the same parameters given in the construction of the simulator
simulator.run()
One simulation has been performed
print(simulator.simulations_performed)
1
The simulation can be accessed at simulator.simulation[0]
simulator.simulation[0].results
{'times': array([ 0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ,
1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1,
2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2,
3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4. , 4.1, 4.2, 4.3,
4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5. , 5.1, 5.2, 5.3, 5.4,
5.5, 5.6, 5.7, 5.8, 5.9, 6. , 6.1, 6.2, 6.3, 6.4, 6.5,
6.6, 6.7, 6.8, 6.9, 7. , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6,
7.7, 7.8, 7.9, 8. , 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7,
8.8, 8.9, 9. , 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8,
9.9, 10. ]),
'x': array([[-1.39320598, -2.01826781, -1.92234229, ..., 0.57260334,
0.47775196, 0.35274937],
[-0.61909625, -0.59499455, -0.63164027, ..., 2.01253144,
0.98817005, 0.05541641],
[-2.51809684, -2.48784868, -3.03441066, ..., 0.11257338,
-0.40405033, 0.1429242 ],
...,
[-0.82107004, -0.51178216, -0.81680132, ..., 0.78387936,
1.1034507 , 1.02430127],
[-0.33465449, -0.46067137, -0.56094613, ..., 0.78900911,
0.32416553, 0.24886806],
[-2.08355921, -2.09995708, -2.18740501, ..., 0.52437001,
0.55692324, -0.24688365]]),
'power': array([[ 0. , 0.27180005, 0.25864001, ..., 0. ,
0. , 0. ],
[ 0. , -0.06782586, -0.05585131, ..., 0. ,
0. , 0. ],
[ 0. , 0.4282953 , 0.66320787, ..., 0. ,
0. , 0. ],
...,
[ 0. , -0.08141434, -0.02097029, ..., 0. ,
0. , 0. ],
[ 0. , -0.0894174 , -0.06826445, ..., 0. ,
0. , 0. ],
[ 0. , 0.29743996, 0.34384373, ..., 0. ,
0. , 0. ]]),
'work': array([[ 0. , 0.01329571, 0.03928595, ..., 0.32540684,
0.32540684, 0.32540684],
[ 0. , -0.00695885, -0.012054 , ..., -0.56530878,
-0.56530878, -0.56530878],
[ 0. , 0.04366937, 0.09401924, ..., 3.97770679,
3.97770679, 3.97770679],
...,
[ 0. , -0.00762188, -0.01437032, ..., 0.66022454,
0.66022454, 0.66022454],
[ 0. , -0.01207446, -0.02011314, ..., -0.54664451,
-0.54664451, -0.54664451],
[ 0. , 0.02338372, 0.05378687, ..., 2.76370786,
2.76370786, 2.76370786]]),
'heat': array([[ 0.00000000e+00, 2.33633892e-01, 1.83280729e-01, ...,
-2.72725620e-01, -2.27688065e-01, -1.54592886e-01],
[ 0.00000000e+00, 4.66032094e-03, -2.59856982e-03, ...,
1.04164682e+00, 5.29106841e-01, 9.75155947e-01],
[ 0.00000000e+00, -2.45718749e-02, 4.92216685e-01, ...,
-4.16009829e+00, -3.56818263e+00, -4.18657183e+00],
...,
[ 0.00000000e+00, 5.08472014e-02, 9.13533239e-03, ...,
-6.44874457e-01, -6.62877499e-01, -6.67933247e-01],
[ 0.00000000e+00, -3.50177862e-02, -5.70389036e-02, ...,
4.58231926e-01, 6.64349461e-01, 7.18072942e-01],
[ 0.00000000e+00, 1.44641835e-02, 7.03300090e-02, ...,
-2.94412106e+00, -2.95907450e+00, -2.27987359e+00]]),
'delta_U': array([[ 0. , 0.24692961, 0.22256668, ..., 0.05268122,
0.09771877, 0.17081395],
[ 0. , -0.00229853, -0.01465257, ..., 0.47633804,
-0.03620194, 0.40984716],
[ 0. , 0.0190975 , 0.58623593, ..., -0.1823915 ,
0.40952416, -0.20886504],
...,
[ 0. , 0.04322532, -0.00523499, ..., 0.01535008,
-0.00265296, -0.00770871],
[ 0. , -0.04709225, -0.07715204, ..., -0.08841258,
0.11770495, 0.17142843],
[ 0. , 0.03784791, 0.12411688, ..., -0.1804132 ,
-0.19536663, 0.48383428]]),
'energy': array([[3.86527358e-02, 2.85582342e-01, 2.61219418e-01, ...,
9.13339519e-02, 1.36371507e-01, 2.09466686e-01],
[3.62719164e-02, 3.39733904e-02, 2.16193462e-02, ...,
5.12609957e-01, 6.99738822e-05, 4.46119079e-01],
[5.76154506e-01, 5.95252003e-01, 1.16239044e+00, ...,
3.93763007e-01, 9.85678669e-01, 3.67289462e-01],
...,
[8.00398293e-03, 5.12293034e-02, 2.76899152e-03, ...,
2.33540660e-02, 5.35102405e-03, 2.95275967e-04],
[1.10671162e-01, 6.35789147e-02, 3.35191174e-02, ...,
2.22585781e-02, 2.28376113e-01, 2.82099594e-01],
[2.93525140e-01, 3.31373047e-01, 4.17642022e-01, ...,
1.13111943e-01, 9.81585096e-02, 7.77359420e-01]])}
Let’s perform another simulation with a different set of parameters.
simulator.run(
tot_sims = 100_00,
dt = 0.01,
tot_steps = 7_000,
snapshot_step = 500
)
Now 2 simulations have been performed
print(simulator.simulations_performed)
2
The second simulation can be accessed at simulator.simulation[1]
simulator.simulation[1].results
{'times': array([ 0., 5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60.,
65., 70.]),
'x': array([[-0.87221322, -1.0955517 , 0.51170871, ..., -0.10011611,
-0.77325558, 0.47818033],
[-2.44897706, 1.22207542, 0.84702523, ..., 1.63530661,
1.34227907, 0.60301154],
[-2.14848429, -0.06627565, -0.2994279 , ..., 0.27821104,
0.14949696, 0.35752 ],
...,
[-0.88761793, 2.94156042, 0.05477785, ..., 1.63768262,
1.1136945 , -0.42792164],
[-0.85540488, -1.44017708, 0.30362297, ..., 2.37718747,
1.35456978, 1.81298381],
[-0.83874356, 1.7709441 , 0.65478371, ..., 2.85621192,
1.18424655, -0.24347578]]),
'power': array([[ 0. , 1.0561501 , 0. , ..., 0. ,
0. , 0. ],
[ 0. , -0.08707466, 0. , ..., 0. ,
0. , 0. ],
[ 0. , 0.48213174, 0. , ..., 0. ,
0. , 0. ],
...,
[ 0. , -0.5881639 , 0. , ..., 0. ,
0. , 0. ],
[ 0. , 1.27201877, 0. , ..., 0. ,
0. , 0. ],
[ 0. , -0.27915072, 0. , ..., 0. ,
0. , 0. ]]),
'work': array([[ 0. , 0.98683144, 0.98683144, ..., 0.98683144,
0.98683144, 0.98683144],
[ 0. , 0.45242369, 0.45242369, ..., 0.45242369,
0.45242369, 0.45242369],
[ 0. , 1.23391027, 1.23391027, ..., 1.23391027,
1.23391027, 1.23391027],
...,
[ 0. , -1.56171708, -1.56171708, ..., -1.56171708,
-1.56171708, -1.56171708],
[ 0. , 1.39418006, 1.39418006, ..., 1.39418006,
1.39418006, 1.39418006],
[ 0. , -0.56298211, -0.56298211, ..., -0.56298211,
-0.56298211, -0.56298211]]),
'heat': array([[ 0. , 1.20475466, -0.87169961, ..., -0.38578608,
0.58130387, -0.85476592],
[ 0. , -0.95264858, -0.96560668, ..., -0.77550008,
-0.91872984, -0.8985074 ],
[ 0. , -0.99519243, -0.71940787, ..., -1.30317466,
-1.2019866 , -1.35727404],
...,
[ 0. , 3.44338808, 2.0052821 , ..., 1.76187921,
1.56502287, 2.57803975],
[ 0. , 1.5778251 , -1.15693651, ..., -0.45108433,
-1.33654713, -1.06893565],
[ 0. , 0.85365861, 0.61606835, ..., 2.27924254,
0.5734546 , 1.32959721]]),
'delta_U': array([[ 0. , 2.19158609, 0.11513183, ..., 0.60104536,
1.5681353 , 0.13206552],
[ 0. , -0.50022488, -0.51318299, ..., -0.32307638,
-0.46630615, -0.44608371],
[ 0. , 0.23871784, 0.5145024 , ..., -0.06926439,
0.03192367, -0.12336377],
...,
[ 0. , 1.881671 , 0.44356503, ..., 0.20016213,
0.00330579, 1.01632267],
[ 0. , 2.97200516, 0.23724355, ..., 0.94309573,
0.05763293, 0.3252444 ],
[ 0. , 0.29067649, 0.05308623, ..., 1.71626043,
0.01047248, 0.7666151 ]]),
'energy': array([[0.00408237, 2.19566846, 0.11921419, ..., 0.60512773, 1.57221767,
0.13614788],
[0.52488363, 0.02465875, 0.01170064, ..., 0.20180725, 0.05857748,
0.07879992],
[0.32975404, 0.56847188, 0.84425644, ..., 0.26048965, 0.36167771,
0.20639027],
...,
[0.00315743, 1.88482843, 0.44672246, ..., 0.20331956, 0.00646322,
1.0194801 ],
[0.00522694, 2.97723209, 0.24247049, ..., 0.94832266, 0.06285987,
0.33047134],
[0.00650091, 0.2971774 , 0.05958714, ..., 1.72276134, 0.01697339,
0.77311601]])}
Let’s analyse the results of the first simulation
sim = simulator.simulation[0]
Plot trajectories number 1, 2, 5 and 10
sim.plot_sim('x', [1, 2, 5, 10])
One can plot the following quantities:
sim.result_labels
['x', 'power', 'work', 'heat', 'delta_U', 'energy']
Plot the work done of trajectories 0 to 9
sim.plot_sim('work', range(10))
Plot the average position, the variance of the position, the average work and heat.
sim.plot_average('x')
sim.plot_variance('x')
sim.plot_average('work')
sim.plot_average('heat')
Plot an animation of the PDF of the position
sim.animate_pdf('x', x_range = [-5, 5], y_range= [0, 0.5])
Plot the PDF of work
sim.animate_pdf('work', x_range = [-2, 7], y_range = [0, 1.5])
To perform simulations with a different potential, a new simulator has to be created. Let us create one for a non harmonic potential.
def double_well_potential(x,t):
"""Quartic potential with a double well with moving center"""
tf = 5.0
v = 0.5
if t < tf:
return (x-v*t)**4 - (x-v*t)**2
if t >= tf:
return (x-v*tf)**4 - (x-v*tf)**2
simulator2 = Simulator(
tot_sims = 100_000,
dt = 0.001,
tot_steps = 10_000,
snapshot_step = 100,
harmonic_potential= False,
potential = double_well_potential
)
simulator2.run()
sim2 = simulator2.simulation[0]
sim2.animate_pdf('x', x_range=[-2, 4.5], y_range = [0, 1])
By default, the initial condition is at thermal equilibrium at t=0. But this can be changed. Let’s perform some simulations with a fixed initial condition.
def init_condition():
"""Initial condition particle is at x=4.0 at time t=0"""
return 4.0
Let’s create a new simulator with that initial condition and the quartic potential
simulator3 = Simulator(
tot_sims = 100_000,
dt = 0.001,
tot_steps = 10_000,
snapshot_step = 100,
harmonic_potential= False,
potential = double_well_potential,
initial_distribution = init_condition
)
simulator3.run()
simulator3.simulation[0].animate_pdf('x', x_range = [-2, 4.5], y_range = [0, 0.8])